
install.packages("foreign"); install.packages("MASS");install.packages("stargazer");install.packages("AER")
library(foreign); library(MASS);library(stargazer);library(AER)




#setwd("~/Dropbox/Job Market Paper/IO Revisions/Replication_Files/Replication_Appendix/A_Fig_13")




dd<-read.dta("data_fig_a13.dta")

ft1<-glm.nb(Entrants~iron,data=dd)

ft1a<-glm.nb(Entrants~iron+X+Y+rivers,data=dd)

ft2<-glm.nb(Entrants~textile,data=dd)

ft2a<-glm.nb(Entrants~textile+X+Y+rivers,data=dd)

ft3<-glm.nb(Entrants~total,data=dd)

ft3a<-glm.nb(Entrants~total+X+Y+rivers,data=dd)




a<-list(ft1,ft1a,ft2,ft2a,ft3,ft3a)


print(
stargazer(list(coeftest(ft1),coeftest(ft1a),coeftest(ft2),coeftest(ft2a),coeftest(ft3),coeftest(ft3a)), keep=c("iron","textile","total"))
)


dispersion<-rep(NA,6)

for(i in 1:6){
	
	dispersion[i]<-a[[i]]$theta
	
	
}



print("This is the Dispersion Parameters")

print(dispersion)


newdata<-data.frame(seq(0,8,by=1))

names(newdata)<-"total"

pdf(file="Figure_A13.pdf")


par(mar=c(3.75,3.75,3.75,1))

plot(newdata$total,predict(ft3,newdata=newdata,type="response"),type="l",xlab="",ylab="",lwd=2)

grid()

names(newdata)<-"textile"


lines(newdata$textile,predict(ft2,newdata=newdata,type="response"),lty=2,lwd=2,type="l")

names(newdata)<-"iron"

lines(newdata$iron,predict(ft1,newdata=newdata,type="response"),lty=3,lwd=2,type="l")


text(6,12,"Total")
text(4,15,"Textiles")

text(5.625,17.5,"Iron")

mtext("Proto-Industrial Centers",1,line=2.2)
mtext("Number of States",2,line=2.2)
mtext("Proto-Industry & The Number of States",3,line=1.5,cex=1.5,font=2)


dev.off()